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The primary motivation for work presented in this paper has been to contribute to the continuing integration 
■ of quantum measurement theory with traditional (engineering) disciplines of measurement and control. Various 
researchers engaged in this endeavor have found that the concepts and methods of theoretical engineering provide a 
fresh perspective on how differences and relationships between quantum and classical metrology can be most cleanly 
' understood. This approach has been especially fruitful in scenarios involving continuous measurement, for which 
a number of important physical insights and results of practical utility follow simply from the formal connections 
between quantum trajectory theory and Kalman filtering |], g, [|, Q |[ [| Q . 

Here we describe a general formalism for parameter estimation via continuous quantum measurement, whose equa- 
tions are amenable to analytic and numerical optimization strategies. In addition to being useful for practical design 
of quantum measurements, we find that this approach sharpens our understanding of the significance and origin of 
I . Standard Quantum Limits (SQL's) in precision metrology. Following the basic notion that the "standard limit" for 
' any measurement scenario should be derivable by optimization over some parametric family of "standard" measure- 
ment strategies, we present results that generalize the SQL for force estimation through continuous monitoring of the 
position of a test mass. Our analysis shows that the canonical expression for the force SQL in continuous position 
measurement stems from a rather arbitrary limitation of the set of allowable measurement strategies to those with 
^ c"| constant sensitivity, and we find that a lower expression (by a factor of 3/4) can be obtained when time variations 
i| are allowed. It follows that further expansions of the optimization space (such as adaptive measurements with real- 
. time feedback ]l|) should be considered in order to arrive at an SQL that consistently accounts for a natural set of 
measurement strategies that are "practically equivalent" in terms of inherent experimental difficulty. 

For clarity, the main results of this paper are presented in the first and third sections within the concrete context 
of force estimation via continuous position measurement. In order to emphasize the general nature of our formalism 
and the conclusions we derive from it, the second section provides a more abstract development that arrives at all the 
equations needed for sensitivity optimization in a broad class of continuous measurement scenarios. As this general 
treatment is rather technical, we note that it is not crucial to the overall logical flow of the paper. Very recently 
Gambetta and Wiseman have discussed a similar approach to parameter estimation for resonance fluorescence of a 
two-level atom paying particular attention to how information about the unknown parameter, and also about the 
quantum state, changes with different kinds of measurements ||. 
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I. FORCE ESTIMATION BY CONTINUOUS MEASUREMENT OF POSITION 



The aim of this section is to present a formalism for continuous parameter estimation in the specific context of 
a harmonic oscillator subject to an unknown force linear in x. This section gives a rigorous and a more general 
treatment of the ideas previously worked out by one of us Q . We first derive the conditional evolution equations for 
the oscillator under continuous position measurement, then discuss their control-theoretic interpretation as Kalman 
filtering equations. We then show how a Bayesian parameter estimator can be obtained from the Kalman filter in 
this scenario. 
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A. Conditional evolution equations 



We will derive the equations of motion of a continuously observed system conditioned on the measurement record. 
Our treatment is based on the model of continuous measurement of Caves and Milburn ||, which in turn was based 
on work of Barchielli et al |l0| . Their derivation is solely based on the standard techniques of operations and effects 
in quantum mechanics which makes it very transparent. Similar results could have been obtained by making use of 
the quantum-stochastic calculus of Hudson as was done by Belavkin and Staszewski jl^] . 

In continuous measurement — often an accurate description of experimentally realizable measurements — projective 
collapse of the wavefunction, and hence also the Zeno effect, can be avoided by continually performing infinitesimally 
weak measurements. A weak measurement consists of weakly coupling the system under interest to a (quantum- 
mechanical) meter, followed by a von Neumann measurement of the meter state. As there was only a weak coupling, 
only very little information about the system of interest is revealed and there will only be a limited amount of back- 
action. At first we will introduce the concept of weak measurements in the framework of position measurement. 
Then we will show how to derive the equations of motion for a quantum particle subject to a whole series of weak 
measurements. The treatment of continuous measurements will then be obtained by taking appropriate limits. 

The aim of a weak position measurement is to get some information out of the system, although without disturbing 
it too much. This can be done by applying a selective POVM {A^(x)} where there is a lot of overlap between the 
A(: (x) associated with different measurement results £. This overlap is proportional to the variance of the measurement 
outcome, but inversely proportional to the variance of the back-action noise. As shown by Braginsky and Khalili [ p^[ , 
the product of those variances always exceeds fi /4. Equality is achieved if and only if A^[x) is Gaussian in x. As 
we are interested in the ultimate limits imposed by quantum mechanics, we will assume our measurement device is 
optimally constructed so as to yield a Gaussian A^(x): 
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This is equivalent to the model of Barchielli and also of Caves and Milburn || who obtained it by explicitly working 
out the case of linear coupling between a (Gaussian) meter and the particle followed by a von Neumann measurement 
on the meter. 

We will now assume that the wavefunction of the observed particle is also Gaussian. This is a reasonable assumption 
as we will soon take the limit of many Gaussian measurements, each of which effects a Gaussian "conditioning" of 
the particle's wavefunction. Ultimately the wavefunction itself will become Gaussian, whatever its original shape. We 
furthermore assume that the Hamiltonian of the unobserved particle would be given by: 

H« = f + ^x 2 + 9x, (1) 
2m 2 

where 9 is the (eventually time-dependent) force to be estimated. It will turn out to be very useful to parameterize the 
Gaussian wavefunction of the particle by a complex mean x = x r + Hi and complex variance a — ay + i&i (throughout 
the paper the notation a instead of a 2 will be used to denote the variance) : 

- ~ &i ~ p. i 

X Xy I ~ X i P ih ~ 

a r <7 r 

H 2 —t -r -j TlGi 



Ax 2 = 7-1- Ap 2 = — - AxAp + ApAx = — (2) 
2oy 2(7,. cr r 

The values of these quantities will in general depend on the value of 9. In this subsection we will supress this 
dependence but in the following we will denote the mean position conditioned on a particular value of 9 by xg and 
likewise for the other expectation values. We will now derive the dynamics of this state if a measurement takes place 
at time r. From time to t", just before the measurement, the equations of motion are governed by the Schrodinger 
equation: 

da ih ( m 2 uj 2 2 \ dx a(t) . 2 , 

1 — T?-v{t) -77 = ~rr \ 9 + mu} x ) ( 3 ) 



dt m V h J dt ih 
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The corresponding x, p and second order moments can easily be derived. The equation for a indicates the spreading 
contracting of the wavepacket induced by the harmonic oscillation. At time r, the POVM {A^(x)} is performed. £ will 
be a Gaussian distributed random variable with expectation value x(t~) and variance D + Ax 2 (t~). Straightforward 
calculations show that the post-measurement wavefunction, conditioned on the result £, is parameterized by: 

x d T ) = ~ /" — \ , n ( 4 ) 



o{t) cr(r-) D 4V ' a{r~) + D 

The equation for a now indicates the contracting effect of the position measurement. The expectation values x and 
p become: 

x(r) = x(r-) + ^f^-(C-x(r-)) 
a r (T)D 

P{r) = p(r-) + -^^(e-S(r-)) (5) 

Note that the back-action manifests itself by constantly introducing white noise, i.e. £ — x(t~), into the system. 

It is trivial to write down the dynamical equations in the case of a finite number (TV) of measurements: we just have 
to repeat the previous two-stage procedure N times. However we are interested in taking the limit of infinitesimal time 
intervals dt between two measurements. This will only make sense if at each infinitesimal time step the wavefunction 
is only subject to an infinitesimal disturbance. Referring to equation this implies that the measurement accuracy 
D has to scale as 1/dt. Therefore we define the finite sensitivity k by the relation D = l/(kdt), implying that only an 
infinitesimal amount of information is obtained in e ach me asurement. In this limit, the random zero- mean variable 
(£ — x(t~))/D has a standard deviation given by ^kdt/2. This is very convenient as a Gaussian random variable 
with zero mean and variance sfdt is by definition a Wiener increment, and therefore we can make use of the theory 
of Ito calculus. Defining d£(t) — ^dt as being the measurement record, and using the notation of Ito calculus, the 
complete equations of motion conditioned on the measurement result for a Gaussian particle subject to continuous 
observation of the position can be written down: 

= x (t)dt + vt(t)dW (6) 

dx(t) = i^-dt + v x (t)dW (7) 
m 

dp(t) = -muj 2 x(t)dt~9{t)dt + v p (t)dW (8) 
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If the sensitivity k is kept constant during the whole observation (Vt, k(t) — k(0)), equation @ can be solved exactly. 
Given initial condition erg, the solution is: 



w 00 y gB±3L exp(2ffif) + 1 y V m V ' 

This shows that the position variance of the wavefunction evolves at least exponentially fast to a steady state. 
The damping is roughly proportional to the square root of the sensitivity, while the steady state solution has a 
variance inversely proportional to it. This result means that a continuously observed particle is localized, although 
not confined, in space. It is interesting to note that this localization increases with the mass of the particle, such 
that it is very difficult to localize a light particle. Indeed the steady state position variance can be understood from 
the point of view of Standard Quantum Limits for position measurement [[l3|. For example if lu 2 3> hk/m then 
Aa; 2 oo — h/2mu>. Similarly, if we take t — l/i?e[f2] to be the time for an effectively complete measurement, then for 
a free particle Ax 2 ^ = ht/m and so the steady state position variance is the same as the SQL for ideal position 
measurements separated by time intervals of length l/i?e[fi]. 
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B. Kalman filtering interpretation 

Let us now try to give a "signal processing" interpretation to equations (p|-[Toj) . The Wiener increment was denned 
as the difference between the actual and the expected measurement result. As it is white noise, it is clear that the 
expected measurement result was actually the best possible guess for the result. This is reminiscent to the innovation 
process in classical control theory: the optimal filtering equations of a classical stochastic process can be obtained by 
imposing that the difference between the actual and expected (i.e. filtered) measurement be white noise. Indeed, in 
a previous paper Q , one of us noticed that the equations (P|jlO|) have exactly the structure of the Kalman filtering 
equations associated with a classical stochastic linear system. This is in complete accordance with the dynamical 
interpretation of quantum mechanics as describing the evolution of our knowledge about the system. 

The classical stochastic system that has exactly the same filtering equations as our continuously observed quantum 
system is given by: 




dVi and (SV2 are two independent Wiener increments and correspond to the process noise and measurement noise 
respectively. It is very enlightening to look at the corresponding weights of these noise processes: the higher the 
sensitivity, the more accurate the measurements, but the more noise is introduced into the system. Moreover measuring 
the position only introduces noise into the momentum. This clearly is a succinct manifestation of the Heisenberg 
uncertainty relation. Indeed, the product of the amplitude of the noise processes of measurement and back-action is 
independent on the sensitivity k and exactly given by fi/2. 

The equations for the means xg and pg are now given by the Kalman filter equations of this classical system, and 
the equations for the variances Axg, Apg, AxgApg + ApgAxg are given by the associated Riccati equations. This is 
very convenient as this will allow us to use the convenient language of classical control theory to solve the estimation 
problem. 

C. Continuous Parameter Estimation 

Let us now consider the basic question of this paper: how can we get the best estimates of the unknown force {0(t)} 
acting on the system, given the measurement record {<i£t}? The natural way to attack this problem is the use of 
Bayes rule. As we have a linear system with {d£t} a linear function of {9(t)}, and the noise in the system is Gaussian, 
this will lead to a Gaussian distribution in {8(t)}. Moreover, due to the linearity, the second order moments of this 
distribution will be independent of the actual measurement record. Therefore the accuracy of our estimates will only 
be a function of the sensitivity chosen during the observation process and of the prior knowledge we have about the 
signal {0(t)} (for example that it is constant) . This will allow us to devise optimal measurement strategies. 

The formalism that we have developed is particularly useful in the case that we parameterize {9(t)} as a linear 
combination of known time-dependent functions {fi(t)}, but with unknown weights {9i}: 

n 

0(t) = X>/i(*) 

i=l 

The estimation, based on Bayes rule, will lead to a joint Gaussian distribution in the parameters {9i}. Indeed, we 
have the relations: 

p({9i}mt+dt)}) ~ p{m\{(>i}dm})p({Oi}\m}) 

~ p(dm\^(tAdi},m}))pm\{m}) (^) 

In the last step we made use of the fact that the Kalman estimate X{g.}(t) is a sufficient statistic for d£(t). Moreover 
all distributions are Gaussian, while x^g.y(t) is some linear function of {9i} due to the linear character of the Kalman 
filter: 



%*}(*)=!> / dt'g(t,t')Mt>) 
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The function g(t,t') can easily be calculated using equations (@|uj|). To obtain the variance of the optimal estimates 
of {9i}, formula ( |l3| ) has to be applied recursively. By explicitly writing out the Gaussian distributions, and making 
use of the fact that the product of Gaussians is still a Gaussian, it is then easy to show that the variances at time r 
are given by: 



1 

<*6i 



dt 

W) 



(14) 



A more intuitive way of obtaining the same optimal estimation, given a fixed measurement strategy, of {9i} can be 
obtained by a little trick: we can enlarge the state vector (xg,pe) with the unknowns, and construct the Kalman 
filter and Riccati equation of the new enlarged system, xo and pg, till now the expected values conditioned on a fixed 
value of the force, then get the meaning of the mean of these expected values over the probability distribution of the 
unknown force. In other words, the new x and p become the ensemble averages over the pure states labeled by a fixed 
force 9. The new enlarged system, in the case of one unknown parameter 9, reads: 




y/2k{t)dVl 



(15) 



(16) 



The Kalman filter equations will give us the best possible estimation of the vector (x,p,9) at each time, while the 
Riccati equation determines the evolution of the covariance matrix P: 



J t I | I = A (t) | P | + 2k(t)P(t)C T | <j£(t) - C | p 
P = A(t)P + PA T {t) -2k(t)PC T CP + 2k{t)BB q 



(17) 



(18) 



An optimal measurement strategy, dependent on the sensitivity, will then be this one that minimizes the (3, 3) 
element in P at time tfi na i- An analytic solution of this problem does not exist in general, as the Riccati equations 
are quadratic. However, in the case of constant f(t) = /(0) and constant sensitivity k(t) — k(0) analytical results will 
be derived. 

Before proceeding however, it is interesting to do a dimensional analysis to see how the variances will scale. We 
begin by scaling t = tjr with r the duration of the complete measurement. Introducing the matrix 



T = 




(19) 



T 1 PT 1 is dimensionless. If we then scale the sensitivity as k(t) = k(t)hT 2 / (2m) , 
= 9 ^Jhm/2T 3 and do the appropriate transformations B — > B and C — > C, we get the equivalent state 



it can easily be checked that P 
the force 
space model: 



.4 



1 

-^T 2 f(t) 





B 



c 



1 



(20) 



The new filter equations are still given by with the substitution (A,B,C,k(t)) -> (A, B, C, k(i)). This 

observation has an immediate consequence if we are measuring the force acting on a free particle (lj = 0): the 
standard deviation on our estimate will always scale like y^hm/2T 3 , and the chosen sensitivity will only affect the 
accuracy by a multiplicative pre-factor. 
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II. GENERAL FORMALISM FOR QUANTUM PARAMETER ESTIMATION 

In this section we develop a description of the problem of estimating unknown parameters 9 of the dynamics of a 
quantum system from the results of generalized measurements. This general problem can be addressed in essentially 
the same way as the specific problem of force estimation for an oscillator that was discussed in the previous section. 
An approach to this problem has been proposed by one of us j| and we will formulate the theory in the language 
of operations and effects and consider in particular the case of measurement currents that are continuous in time, as 
in the case of homodyne detection pi or continuous position measurement. The fundamental basis of this approach 
is to propagate an a posteriori probability distribution p (#|I[ 0t )) for the parameter 9 conditioned on the history of 
measurement results I[o,t) up to time t by employing Bayes' rule and using the theory of operations and effects to 
calculate the relative likelihood of the known measurement record as a function of 9. Readers who are less interested 
in mathematical details and more interested in the application of our formalism to the force estimation problem may 
skip this section. 



A. General Theory 



We will treat the quantum parameter estimation as an essentially classical parameter estimation problem coupled 
to the quantum measurement updating rules. For each value 9' of 9 there will be a conditioned state pg> describing 
the state of the quantum system conditioned on the measurement history and a particular value of the unknown 
parameter 9. This density matrix would be our best description of the state if we knew the measurement record and 
also that 9 took this particular value. However the value of 9 is not assumed to be known exactly and is described by a 
probability distribution p(9). Hence the density matrix describing the state from the point of view of the experimenter 



p = J d9p(9) Pe . (21) 



The most general quantum evolution and measurement can be described by the theory of operations and effects. 
The following discussion will adapt the treatment of Wiseman and Diosi to our problem . In this work we assume 
that either the dynamics or the measurement are unknown and belong to a family parameterized by 9. Thus we 
consider quantum measurements characterized by a set of operators f2e jT . where 9 labels the value of the unknown 
parameter and r labels the measurement result. Thus there is a separate measurement for each value of 9 and the 
operators are constrained by completeness 

dAt0,o(r)l4 >r Oe, r = 1. (22) 

Here dpepir) is a normalized measure on the space of measurement results r. As in the standard theory, the probability 
of the measurement result r conditioned on 9 is 

dp e (r) = dp e ,o( r ) Tr ^\^e, r pe ■ (23) 

The state of the quantum system after the measurement conditioned on the pair {9, r) is 

dfie t0 (r)Q,e, r petig r fl e , r pe^l r 



dp g (r) 



Tr 



(24) 



If the result of the measurement is unknown or disregarded then the state of the system is an average over the 
conditioned states weighted by their probabilities 

Pe= dp, g (r)p' gr = \ dfj,e i0 (r){le,rpefig ir - (25) 



This is the state of the system conditioned on a particular value of 9 but not on any measurement result. 

The unconditioned probability of the measurement results is found by averaging over the probability distribution 
for 9 and is given by the measure 

dp(r)= [ d9p{9)dp e {r) = [ dp(9)dp g (r). (26) 

J J 9 
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After the measurement we will require that the state conditioned on the measurement result r but not on the value 
of 9 may still be written in the form of Eqn. ( |2l| ) as an average over the states conditioned on particular values of 9, 
thus 

' dp r (d)p' e>r (27) 

for some measure dp r (9) on the space of possible 9. This new measure describes the probability of 9 conditioned on 
the measured value of r. This conditioned probability distribution for 9 is precisely what we wish to calculate. For 
consistency it must be the case that if the measurement result is unknown or disregarded the appropriate state is 
again an average over the conditioned states 

P = / dfx(r)p' r = / dp{r)dp r {9)p' g (28) 



In order to calculate dp r {9) we need to develop a Bayes rule that relates all the probability measures we have 
introduced. In order to do this we note that p' must also be able to be expressed as an average over the probability 
for 9 of the states p' a , thus 



p' = / dp{9)p' e = / dp(9)dp e {r)p' s . r . (29) 



This leads us to the Bayes rule 



dp(r)dp r (9) = dp(9)dp e (r) (30) 



which allows us to calculate dp r {9) in terms of dp{9), the measure that characterizes our prior knowledge about 9, 
and the measures dpe(r), which are part of our specification of the parameterized family of measurements. 

In principle this allows us to optimally update the probability distribution for the unknown parameter in any 
quantum measurement. We are most interested here in the case of measurements that are continuous in time. In 
this situation we wish to derive a stochastic differential equation that updates the distribution for 9 conditioned on 
measurement current. Since the case of photon detection style measurements is considered in |}| we will consider 
measurements like homodyne detection where the measurement results are continuous but not differentiable functions 
of time, in |l5[ ] these are termed diffusive measurements. This will require that we develop stochastic differential 
equations to describe the measurement process. 

For simplicity we will consider the case where there is only a single measurement being made and we will describe 
the measurement result r in an infinitesimal time interval [t,t + dt) by the complex number I(t). We define the 
measurement operators 

n B ,i = l-iH{e)dt - ~tfc + I*cdt. (31) 

These measurement operators may be derived, for example, as the continuous limit of a model of repeated measure- 
ments [^j or from models of quantum optical measurements such as heterodyne or homodyne detection [ jHJ . For 
simplicity we consider the case where there is only one measurement current, the general case may easily be treated 
following the formalism of We also assume that the specific measurement that is being made is known (that is 
that the operator coupling the system to the bath and the measurement made on the bath are known) and so 9 only 
parameterizes the Hamiltonian evolution of the system. This is the most interesting case and simplifies the treatment. 
The extension to the case where the measurement is known but the free system evolution is not unitary but is rather 
described by a Markovian master equation is also straightforward. Now the measurement operator is constrained by 
the completeness relation Eq. (p2|) and this requires that 



dpL efi {I){Idt) = (32) 

dp gfl {I)(I*dt){Idt) = dt. (33) 

These moments mean that we may identify Idt as a complex Wiener increment under the measure dpgfi(I)- However 
in order to specify this measure completely we must also specify the remaining second order moment of the Wiener 
increment (clearly this must also be of order dt) . We will say that 

dp 9fi {I)(Idt)(Idt) = udt, (34) 
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where we need \u\ < 1. In line with our assumption that the measurement interaction and the measurement on the 
bath is known we will require that u is independent of 9. The case u = corresponds in the quantum optical setting 
to heterodyne detection, while \u\ = 1 corresponds to homodyne detection with some local oscillator phase. Note that 
these moments are independent of 9 and so we can drop the subscript 9 for this measure on / from here on. Since the 
moments of Idt under dpo(I) indicate that we consider Idt to be a complex Wiener increment, we adopt the Ito rules 



(Idt) 2 = udt, (Pdt)(Idt) = dt. 



(35) 



Now we would like to know the observed statistics of I under the physical measure dp(I). There are two kinds of 
conditioned expectation values for operators a in this problem. Expectation values conditioned on a particular value 
of the unknown parameter will be denoted dg —Tv[dpg]. On the other hand expectation values conditioned only on 
the history of measurement results will be denoted a =Tr[dp\. Now we know from the preceding discussion that 



dp(9)dp (I)Tr SV e jSl e jp e 

dpo(I) I dp(9)Tr [(1 + I* cdt + IcUt)f 
Je 

dp (I)(l+I*dtc + Idtc i ). 



Hence the expected value of I is 



{I) 



dp(I)I — u<y + c. 



(36) 

(37) 
(38) 

(39) 



From Eq. ( p8| ) we can see that the second order moments of Idt are independent of the state and of 9 and are 
equal to the second order moments under dpo(I)- Thus the transformation from the measure dpo(I) to dp(I) is a 
transformation of drift similar to a Girsanov transformation J16| and we can identify Idt with 



Idt = uc f + cdt + dW 



(40) 



where dW is a complex Wiener increment under the measure dp(I) obeying dW 2 = udt, dW*dW = dt. 

On the other hand the probability measure for the measurement trajectories conditioned on a given value of 9 is 



dp. e (I) = d/J, (I)Tr tigjSlgjpg 

= dm>(I){l + I* dtc e + Idtcl) 



(41) 
(42) 



Using Eq. ( |30| ) it is now straightforward (keeping terms up to second order in Idt) to update the probability for 9 
conditioned on / 

d Mi [Q , f+ dt) = 1 + (eg - c) (l*dt - u*cdt - c ] dt) + (d - c f ) (Idt - cdt - uc^dt) dpi [o t) {9). (43) 

This allows us to write down a stochastic Fokker-Planck equation for the probability distribution of 9 

dp(9)\I [0 , t+dt) ) = [(eg - c) (Pdt - u*cdt - cUt) + H.c] p(0|I [o , t) ) (44) 

Note that under dp(I) the innovation Idt — cdt — uc^dt is a Wiener increment and thus has mean zero and is not 
correlated with either the quantum state or p(9). This equation is very similar in form to the Kushner-Stratonovich 
equation that arises in classical state estimation problems [|l8| . In order to be able to propagate this equation for the 
probability distribution of 9 we must also be able to update the conditioned state pg and hence the expectation values 
eg. From Eq. (pil) we can show that pg obeys the stochastic master equation (SME) 



dpg 



-i[H (9) , pgjdt + V[c]p e dt + H[c (Pdt - cUt - u*c e dt)] pg. 



(45) 



Equation (Q) and the family of stochastic master equations ( [45]) describe the quantum parameter estimation 
problem for measurements with continuous measurement currents such as optical homodyne detection. As we indicated 
at the start of this section, and as in the algorithm discussed in B, a family of quantum states conditioned on 
the measurement record and on different values of 9 is propagated using appropriate SME's while the conditioned 
probability distribution for 9 is propagated using a stochastic Fokker-Planck equation of the kind that arises in 
classical estimation problems. As we shall see below it is possible to solve these equations for certain linear models 
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such as force estimation due to position measurement on a free particle or oscillator. In general it will be necessary 
to integrate these equations numerically after first discretizing 9. In principle this is straightforward although the 
discretization must be sufficiently fine that a good approximation for the mean c + uc' is maintained at all times and 
this will usually involve a prohibitive computational cost. One way of avoiding this is to consider a linear variant of 
this update equation which is in fact more closely allied to the algorithm in This variant is an analogue both 
of the linear version of the stochastic master equation []l7] and of the Zakai equation which is the linear counterpart 
to the Kushner-Stratonovich equation Jl8| | in classical state estimation. This linear variant does not preserve the 
normalization of p (8\I) but does not depend on uc + and yet still propagates the relative probabilities of different 
values of 9. 

The basic observation is that in the Bayes' rule Eq. ( |30| ) the measure d/J,(r) is independent of 9 and only ensures the 
normalization of dp r {9). If we are only interested in the relative likelihood of different values of 9 we may consider 
unnormalized measures dp r {9) on the space of possible 9 and replace dp(r) by any measure on r independent of 9. In 
particular for our example of continuous measurements we may choose 

dpi [0:t+dt) (dWo(I) = dMl)#I[o,t)(0)- (46) 

Substituting from Eq. ([l||) we get 

dp (9\l [0 , t+dt) ) = {c e Pdt + c\ldt) p (9\I [0 , t+dt) ) ■ (47) 

Under this linear propagation equation the dynamics of the unnormalized distribution p (6\lm t -\j may be calculated 
for each value of 9 independently. This will make it possible to calculate relative probabilities of a discrete set of 
possible values of 9 given a particular sequence of measurement results with no constraints on the discretization of 9. 

This formalism for the estimation of a classical parameter in quantum dynamics may readily be generalized to 
the case where there is more than one unknown parameter or where the parameter undergoes some known time 
dependence as in the previous section. Another interesting situation that may be treated straightforwardly in this 
formalism is correlating the measurement results from two quantum measurements both of which depend on 9. Here 
we have assumed that apart from the measurement the dynamics of the quantum system is unitary. If this is not true 
(as is the case for less than perfectly efficient detection for example) then it is straightforward to show that the first 
term of Eq. ( fl5| ) is simply replaced by a Liouvillian term describing the noisy dynamics of the system, thus 

dpg = C(9)p g dt + V[c]p e dt +H\c (l*dt - c\dt - u*c g dt) } pg. (48) 



In the next section we will return the problem of force estimation through continuous position measurement of an 
oscillator. We will be most interested in finding the optimum (possibly time-dependent) sensitivity of the measurement. 



B. Force Estimation through Continuous Position Measurement 

The general formalism of this section may be reduced to the parameter estimation problem we considered at the 
start of the paper in the important case of force estimation though continuous position measurement of an oscillator 
(c = \/2kx, u = 1, H{9) = p 2 /2m + mu> 2 x 2 /2 + 9x). In this case it is possible to solve the system of equations ( |45| ) 
and (Q) explicitly. We have the system of equations 

dpg = -i[p 2 /2m + mu 2 x 2 /2 + 9x,p e ]dt + 2kV[x}pgdt + V2kH [x] pg (idt - 2V2kx e dt S j (49) 

dp{9\I [0 , t +dt)) = 2V2l(xe-x)(jdt-2V2kxdtjp(6\I [0tt+dt) ). (50) 

This linear system preserves Gaussian quantum states of the oscillator and Gaussian probability distributions for 9. 
As a result we only need to find stochastic equations for the first and second order moments of the pg and P(9\I). The 
procedure is to apply standard master equation techniques combined with the Ito rules for stochastic differential 
equations to find equations for the moments of x and p, conditioned on a particular value of 0, from Eq. (f49|) as was 
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done in 0. The unconditioned moments result from averaging over p(6\I) 

Ax 2 = J d9p{8)Tr (x - x g f p g 

AxAp = / d9p(6) (Tr [(xp + px)p g ] /2 - x e p e ) 



Ap 



2 - 



AxAd 



ApA6 



d9p(9)Tr (p-pe) pe 
d9p(9)9x g ) - 9x 



A9 2 



d9p{9)9p g j - 9p 
d9p(9) (9~9) 2 . 



(51) 
(52) 
(53) 
(54) 
(55) 
(56) 



These moments form the covariance matrix P and it is a straightforward though tedious exercise to show that it obeys 
the matrix Riccati equation (|l8| ) we derived in the first section. Similarly the first order moments obey the equations 
( |l7| ) of the Kalman estimator. 



III. STANDARD QUANTUM LIMITS 

The preceding sections dealt with the problem of optimal estimation of parameters of the Hamiltonian given a 
system that is continuously observed. In this section we will derive the explicit equations of the variances on these 
estimates. 



A. Detection of stationary signals 

Let us first introduce the idea of the standard quantum limit in the context of von Neumann measurements. The 
idea is that a particle is prepared in some optimal way at time 0, such that at time r a projective measurement is 
performed to determine the displacement associated with the force. The optimal preparation is crucial as it has to 
balance the position and the momentum uncertainty. The optimal preparation leads to the expression of the Standard 
Quantum Limit. Consider a free particle with a Gaussian wavcfunction (x\ip) and initial parameters x(0),a(0) (see 
equation (||)) and subject to an unknown force 9. The integrated equations of motion (^) are given by: 

x(t) = xq + 9 (tcr(O)M + t 2 /2m) a(t) = a(0) + i—t 

m 

Suppose that at time r we perform a von Neumann measurement of the position. The probability distribution 
associated with this measurement is given by: 

i 2 " 

2m , 



/ I rr. ft" 



p{x\9) ~ exp 



V 



(57) 



Using Bayes rule assuming a flat prior distribution for 9 the variance on the estimate of 9 given the measurement 
result x can easily be derived: 

2m 2 |qffl| 2 2m 2 (q 2 (0) + (q t (0) + |f) 2 ) 

aB - a r (t)t± ^W 4 ( } 

This function is heavily dependent on the initial conditions of the wavefunction of the particle. The standard quantum 
limit can now be derived by choosing the initial conditions such that a g is minimized. This variance can in principle 
go to zero if we allow (AxAp) to be negative, but we will not consider such "contractive" states |2^, ^l| here. We 
therefore impose the condition Ui (0) > in order to focus our attention on the specific issue of sensitivity optimization. 
The optimal i?(0) is then given by <r(0) = ht/m, and this leads to the expression of the Standard Quantum Limit: 

o-e = (59) 
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It is clear that the square of the amplitude of a detectable force has to be bigger than the variance on its estimation 
to be detectable. Therefore the previous formula is the expression of the minimal force that can be detected by a free 
particle of mass m over a time t. Note that the derived formula exceeds the normal equation of the SQL by a factor 
8 as the standard equation is not derived in the context of parameter estimation. 

We will now apply an analogous reasoning to a quantum particle subject to continuous measurement. The explicit 
expression of the variance on the estimated force was given by equation ([H]). As noted at the end of the first section, 
the resulting variance will be given by the standard quantum limit multiplied by a certain factor. From here on we 
will therefore work in the dimensionless picture as defined in (|20|) . In general it is very hard to find the explicit 
expression for the autocorrelation function g(t,t') in equation (|14|). Things get much more feasible if we do not vary 
the sensitivity during the measurement as the system then becomes stationary. It follows that we can assume that 
the values of the variances reached their steady state values given by equation (|TT|) . After some straightforward linear 
algebra, the explicit expression for g(t,t') in the case of steady state is given by: 

g(t,t') = -rexp(-a(t-t'))sm(b(t-t')) (60) 



6 




LOT 



\ 




(61) 



(62) 



Due to the stationarity of the variances, the autocorrelation function g(t, t 1 ) is indeed only dependent on (t — t'), and 
from here on we will therefore use the notation g(t, t') = g(t — t'). The full expression of the variance on our estimate 
now becomes: 

y = 2kj^ dt(J dt'g(t - t')f(t')\ (63) 

The force that acted on the system was assumed to be of the form 8(t) — 8f(t) with f(t) a known function. Note 
that this expression is dimensionless and has to be multiplied by We next introduce F(ui) and G(ui) the Fourier 
transforms of the functions f(t) ■ um t -n(t) and g(t) ■ uro,i](i), where uro,i](i) is the window function over the interval 
[0, 1]. The damping effect due to the back-action noise is responsible for broadening the spectrum of the harmonic 
oscillator with a width of approximately k/ {lot). Basic properties of Fourier transformations lead to the expression: 

1 1 d^expf^-^) ^ G(ui)G* (lu 2 )F(lji)F*(lo2) (64) 



This formula clearly shows that only the frequencies of the signal F((3) near to the natural frequencies of the oscillator 
G(/3) will be detectable. 

Now we shall explicitly calculate the value of <jg in some different cases. Let us first of all assume that the spectrum 
F(j3) is almost constant for all values where G{j3) is substantially different from 0, i.e. around (3 ~ (<-ot). This is 
realistic in some scenarios of interest for the detection of gravitational waves E3|. Let us furthermore assume that 
lot 3> 1, which means that the period of the oscillator is much smaller then the observation time. Next we observe that 
we are allowed to approximate the sinc((wi — lo2)/2) function by a delta-Dirac function if the width of the spectrum 
G((3), determined by the number fc/(«r), is much bigger then 1. This leads to the expression: 



1 k\F(LOT)\ 2 



— OO 



fc|F(cr)| 2 f°° ^ I 



dco\G(Lo)\ 2 (65) 



d^ TZd I u9 ( 66 ) 



2tt J_ oq (a 2 +b 2 - lo 2 ) 2 + Aa 2 b 2 

\F(b)\ 2 ( 2k/(LOT 



X -^=^=L== (67) 



4:Lot \ ^1 + (2k/ (lot) 2 ) 



12 



The function introduced in the last line is only dependent on 2fc/(ojr) 2 , which can be tuned freely by changing the 
value of our sensitivity. The function x( x ) reaches its maximum value 1 for small values of x, meaning that optimal 
detection requires k <§; (ur) 2 . The derivation however required that 1 <C k/(cjr). Therefore, the optimal choice of 
the sensitivity will be given by a value (wt) <C k <C (lot) 2 , leading to the variance on the estimate: 

9 Alot hm 1 2hmuj 

5- ( 69 ) 



~ \F(b)\ 2 2r 3 \F(b)\ 2 t 2 

This corresponds exactly to the expression of the standard quantum limit for an oscillator []l3f . A similar expression 
can be obtained by explicitly integrating ([33]) with f(t) — S(t). The conditions under which this SQL can be reached 
are: 1. The total duration of the measurement is much bigger then the period of the oscillator; 2. The spectrum of 
the signal to be detected is flat around the natural frequencies of the observed oscillator. 

We will now investigate what happens if this second condition is not fulfilled. In the extreme case, the force to be 
detected is constant, corresponding to a delta-Dirac function in the frequency domain. Again under the condition 
that 1 <C lot -C k/(uir), a good approximation of equation ( |64| ) becomes: 

i.^r^ 2fc /(-) 2 , ( 70 ) 

os {lot) 2 1 + (2fc/(wr) 2 ) 2 

The optimal sensitivity is now given by 2k = (lot) 2 , indicating that one has to choose a much higher sensitivity to 
detect constant forces then resonant oscillating forces. The expression for the SQL for detecting constant forces with 
a harmonic oscillator therefore becomes: 



. 9 hm mhuj 

<J8~2(lot 2 — = 71 

2fi T 



It is now natural to look what happens in the limit of to — > 0, it is if the observed particle is free and only subject 
to a constant force. In that case the explicit integration of ( |62|) becomes possible, as a and b both become equal to 
the sensitivity \fk. Straightforward but long integrations lead to: 

8fc 3/2 

°e = 7 V (72) 

AVk -5 + 8 exp(-Vfc) cos(Vk) - exp(-2Vfc) (2 + cos(2y/k) + sin^) 

Minimization over the sensitivity leads to an expression for the SQL for the detection of a constant force with a free 
particle subject to continuous observation: 

ag ~ 3 — ;-. (73) 

T 6 

Note that this expression differs from the corresponding one derived in [Q, where calculations were done without 
properly accounting for the damping effect of measurement back-action. Comparing this result with (|5^) , the variance 
of our estimate obtained by continuous measurement is 3 times bigger then if we were doing projective measurements. 
This is caused by two factors: at the end of the continuous measurement, there is still a lot of information encoded 
about the force in the wavefunction as the variance on the position at time r is not at all equal to oo. Secondly, the 
previous result was obtained by assuming that the variances of our Gaussian wavefunction were in steady state, and 
this is not necessarily the optimal initial condition. Indeed, it turns out that the optimal initial state (not considering 
contractive states) of the continuously observed particle is a Gaussian state with well defined momentum ((Ap 2 ) <C 1) 
and therefore undefined position (Ax 2 ) ^> 1. This makes sense as the force to be detected can only be seen because 
it manifests itself through the momentum. The fact that the position uncertainty is very large is not so bad as the 
position is continuously observed such that it becomes well-defined very quickly. The expression for the variance on 
the force estimate using this optimally prepared initial state can now be calculated exactly by explicitly solving the 
Riccati equations (|is|): 

2k 3 / 2 (sinh(2Vl) + sin(2\/ft)) 

00 ~ jfe(sinh(2v / jfe) + sin(27fc)) - (cosh(2v / A?) - cos(2a/A;)) 

Optimization over the sensitivity leads to an enhancement of 2/3 in comparison with the steady state case. An even 
bigger gain would have been obtained if a projective measurement at the end of the continuous observation were 
allowed. A realistic way to implement this would be to make the sensitivity very large at the end of the measurement. 



13 



If the matrix -P(l) is the solution of the Riccati equation (|l^) at time t = 1, some straightforward calculations show 
that a projective position measurement reduces the estimator variance by P 2 3 ^/-P^i). The optimal initial conditions 

are still given by ((Ap 2 ) <C 1) and (Ax 2 ) ^> 1. The exact expression of the variance on the estimate in function of 
the sensitivity k is then given by: 

4fc 3 / 4 (cosh(2Vfc) + cos(2%/fc)) 

° e ~ k{cosh(2y/k) + cos(2\/F)) - (sinh(2\/fc) + sin^v^)) 

Minimization over the sensitivity leads to the equation: 

a e ~ 0.752^ (76) 

Therefore we have modestly beaten the usual standard quantum limit by optimally preparing the Gaussian wavepacket 
and doing a von Neumann measurement at the end of the continuous measurement. This shows that a continuous 
measurement together with a projective measurement at the end on a optimally prepared state can reveal more 
information than only projective measurements. In other words, the balance information gain versus disturbance is 
a little bit in favor of continuous measurement. Although noise is continuously fed into the system by the sensor, we 
can extract more information about the classical force. 

An even better performance can be obtained if we vary the sensitivity continuously during the measurement {sen- 
sitivity scheduling) . It is indeed the case that backaction noise introduced in the beginning of the measurement does 
more harm than backaction noise at the end of the measurement, as the random momentum kicks delivered at any 
given time corrupt all subsequent position readouts. In terms of systems theory, the optimal sensitivity as a function 
of time can be found by solving an optimal control problem associated with equation ( |18|) . This optimal control can be 
determined by solving a Bellman equation by using techniques of dynamic programming p!q ] . Due to the nonlinearity 
of the Riccati equation however, this cannot be done analytically. The optimal sensitivity at time t however can 
easily be obtained: it tends to a delta-Dirac impulse such as to mimic a projective position measurement, reducing 
the variance with P 2 3 Defining the cost-function K = P(3,3)(t) — P? 3 1 \(t)/P(i i i)(t), the optimal control 

problem is then well defined and can be solved numerically. Therefore we discretize the total time in for example 
50 intervals, and in each interval we assume the sensitivity has a constant value kj. The solution can then found by 
applying some kind of steepest descent algorithm over these 50 variables {kj}. It turns out that the optimal k(t) 
in the case of a free particle (u = 0) is a smooth monotonously but slowly increasing function of time. In this free 
particle case, the optimal time-varying sensitivity only leads to a marginal gain: the numerical optimization shows 
that the variance of the estimate becomes exactly equal to a factor 3/4 of the usual Standard Quantum Limit (|E^) . 
Nevertheless, we can present this result as a generalization of the usual SQL to include strategies with sensitivity 
scheduling: 

*» = 3J P W 

Much greater improvements can be expected from the application of sensitivity scheduling to the case of a continu- 
ously observed harmonic oscillator. Indeed, the variance on the position of such a particle is small in the middle of the 
well and at the borders, while it is big elsewhere. Therefore the sensitivity should be varied in a sinusoidal manner, 
such as to measure more precisely at the positions where the variance is small. The optimal variation of sensitivity in 
time could be determined by solving a similar optimal control problem to the one explained in the previous paragraph. 
In the limit where projective measurements are allowed, one expects that the optimal variation of sensitivity should 
correspond to stroboscopic measurement ]13J, which is indeed well-known to beat the usual standard quantum limit. 



B. Detection of non-stationary signals 



The techniques introduced in the previous sections can also be used for the estimation of non-stationary signals, 
as one would have for example in the problem of gravitational wave detection when the arrival time of the signal is 
unknown. Suppose for example that we know that the signal to detect is of the form 9{t — t\) = 9of{t — t\) with f(r) 
known but amplitude 9q and arrival time t\ unknown. An effective non-stationary measurement strategy can in fact 
be implemented by constructing a Kalman filter for system ([l2]) assuming that 9 = (assuming /(t) = for t < 0). 
At times t < ti, the quantity d£ — x(t)dt is by construction white noise with variance dt/2k(t). From time t > t\ on 
however, the force will bias this white noise by an amount J t dt'g(t,t')9(t') as the 9 = Kalman filter models the 
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wrong system. This bias will be detectable once it transcends the white noise at time t\ + At: 



ti+At r t I ftt+At ^ 



dt^dt'g(t,t')6 f(t'-h)>^^ ^) CM 

The goal is now to make this At as small as possible. The previous equation can again be solved analytically if one 
has a constant sensitivity and steady state conditions. To make things easier we assume that the observed particle 
is free {ui = 0), although all calculations can be performed in the more general case too. Let us first assume that 
the signal to detect is a kick at time t\\ f(t — ti) ~ 5(t — ti)r with r some measure of the duration of the kick |fl3| . 
Introducing the dimcnsionless parameter k — Aty /, frk/2m, the previous inequality becomes: 

1 Tim k , . 

#0 > -\\— K 7 77 r~s (79) 

t V At 1 — exp(— k)(cos(k) + sm(Kj) 




(80) 

In the last step the optimal k, related to the optimal sensitivity k, was chosen. The meaning of this equation is clear: 
a kick with an amplitude 9q will only be observed after a time span At = 4hm/T 2 0Q. Moreover, the sensitivity has to 
scale inversely with the square root of At. 

An analogous treatment applies to the case of a constant force f(t — ti) = W[o,oo](^ — ^i)- I n this case inequality 
((78|) becomes: 



" V At 3 exp(-K) cos(k) + k - 1 ( ' 

>- < 82 > 

As expected, we recover the well known standard quantum limit, but now in a different set-up. 

The previous arguments can be refined by using techniques of classical detection theory such as the concept of the 
matched filter. The results will however be qualitatively similar to the previous ones. 

More advanced detection schemes can also be constructed by adoptively changing the sensitivity as a real-time 
function of the measurement record Q . A possible application of this is a scheme for the detection of a signal with 
unknown arrival time: first one chooses the optimal sensitivity for estimating the arrival time, and from the moment 
on the signal is detected the sensitivity is brought to its optimal value for detecting the amplitude of the signal. More 
sophisticated versions of this adaptive measurement could be very useful in realistic stroboscopic measurements where 
the initial phase of the harmonic oscillator is unknown, as the measurement sensitivity could be made a real-time 
function of the estimated particle position. 
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